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ABSTRACT 

We present the results of standard one- dimensional test problems in relativis- 
tic hydrodynamics using Glimm's (random choice) method, and compare them 
to results obtained using finite differencing methods. For problems containing 
profiles with sharp edges, such as shocks, we find Glimm's method yields global 
errors ~ 1 — 3 orders of magnitude smaller than the traditional techniques. The 
strongest differences are seen for problems in which a shear field is superposed. 
For smooth fiows, Glimm's method is inferior to standard methods. The location 
of specific features can be off by up to two grid points with respect to an exact 
solution in Glimm's method, and furthermore curved states are not modeled 
optimally since the method idealizes solutions as being composed of piecewise 
constant states. Thus although Glimm's method is superior at correctly resolv- 
ing sharp features, especially in the presence of shear, for realistic applications in 
which one typically finds smooth fiows plus strong gradients or discontinuities, 
standard FD methods yield smaller global errors. Glimm's method may prove 
useful in certain apphcations such as GRB afterglow shock propagation into a 
uniform medium. 

Subject headings: hydrodynamics — methods: numerical — relativity 



1. Introduction 

Interest in relativistic hydrodynamics has heightened in recent years due to the explosion 
in the field of gamma-ray bursts (GRBs — Costa et al. 1997, van Paradijs et al. 1997, Frail 
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et al. 1997, MacFadyen & Woosley 1999, Aloy et al. 2000, Frail et al. 2001, Fox et al. 2005, 
Gehrels et al. 2005, Bloom et al. 2006, O'Brien et al. 2006). The current paradigm for GRBs 
involves the extraction of energy from a newly formed ~ IOMq black hole and coUimation 
into a relativistic jet, which then propagates along the line of sight to the observer. The 
emission is thus strongly beamed and Doppler boosted. The interaction of the jet with 
the circumstellar medium produces afterglow. For Newtonian hydrodynamics the density 
contrast across a strong shock is given by Pshock/Pbackground = (r + l)/(r — 1), where F is 
the polytropic index; in relativistic hydrodynamics Pshock/Pbackground = + l)/(r — 1), 
where 7 is the Lorentz factor (Blandford & McKee 1976). For putative values 7 ~ 10^ — 10^ 
thought to be required for GRB jets, a relativistic shock can have extremely high density 
and be very narrow due to Lorentz contraction. This poses a severe test for standard finite 
difference (FD) methods, and necessitates adaptive mesh refinement (Zhang & MacFadyen 
2006, Morsony, Lazzati, & Begelman 2007). Adaptive refinement techniques also present 
challenges, as it has yet to be demonstrated that increased levels of refinement on a complex, 
multidimensional problem, lead to convergent solutions. The desired test of showing that 
a standard performance metric integrated over the computational volume asymptotes to a 
constant value with increasing level of refinement has yet to be carried out (e.g., Zhang & 
MacFadyen 2006). 

Traditional methods for calculating hydrodynamical evolution of a relativistic fluid have 
relied on finite differencing, i.e., discretizing the differential equations (Norman Sz Winkler 
1986, Marti & Miiller 2003, Del Zanna & Bucciantini 2002, Lucas-Serrano et al. 2004). 
Figure [T] presents an example of smearing inherent in standard FD methods. It shows the 
evolution of Lorentz factor 7 in a spherical relativistic blast wave calculation initialized with 
a Blandford & McKee (1976) solution, taking 70 = 5 initially. Each panel shows the same 
initial conditions, with increasing grid resolution along the +0;— direction. We use a three 
dimensional Cartesian grid and utilize the method described in del Zanna & Bucciantini 
(2002). Our implementation of their method is detailed in Cannizzo, Gehrels, & Vishniac 
(2004). Within each panel the number of grid points along the direction of propagation 
is increased by a factor of 4. In the fourth panel, for which there are 64 grid points per 
small tick mark, one can see the clear development of a forward/reverse shock feature. The 
inherent smearing behavior of the technique is evident by comparing successive panels. 

2. Background 

Glimm (1965) presented the theoretical basis for the random choice, or Glimm's method. 
It relies on first idealizing the solution in (P, p, v) over N grid points as consisting of N 
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piecewise constant states, and then solving the local Riemann problem N — 1 times between 
adjacent grid points. Second, a random location is selected within a cell, the exact solution 
evaluated at that point, and then that is used as the starting solution for the next time step. 
Chorin (1976, 1977) developed Glimm's method into a numerical algorithm for problems 
that could be formulated in terms of nonlinear hyperbolic conservation laws. Sod (1978) 
reviewed several techniques for Newtonian hydrodynamics and found Glimm's method to 
be superior in terms of preserving the sharpness of shock edges. In the early studies using 
Glimm's method one sees clear deficiencies in the solutions, however, both in terms of shock 
front localization and overall stability. 

A breakthrough came from Colella (1982) who proposed using the van der Corput 
sequence instead of a standard random number generator for determining the solution eval- 
uation location with cells in each time step. This sequence is generated by a simple manip- 
ulation of the digits in the binary representation of consecutive integers. 

The application of Ghmm's method to relativistic hydrodynamics became possible when 
Balsara (1994) and Marti & Miiller (1996) generalized the solution of the Riemann solution 
for relativistic hydrodynamics. Wen, Panaitescu, & Laguna (1997) used the results of Marti 
& Miiller to implement a relativistic hydrodynamics Glimm's method. Their study and 
Panaitescu ct al. (1997) are the only works to date that employ Glimm's method for rela- 
tivistic hydrodynamics. 

3. Methodology 

The basic method is demonstrated in Wen et al. (1997, see their Fig. 2). In one half time 
step the exact solution to the local Riemann problem is calculated between two grid points, 
at a position determined by the van der Corput sequence. As explained in Colella (1982), the 
sequence is determined by taking the binary representation of the positive integers, 1 = I2, 
2 = IO2, 3 = II2, 4 = IOO2, 5 = IOI2, 6 = IIO2, 7 = III2, 8 = IOOO2, etc., and then 
flipping the binary digits with respect to the (binary) decimal point, yielding the sequence 
Qi = .I2 = 0.5, 02 = .OI2 = 0.25, as = Ah = 0.75, = .OOI2 = 0.125, 05 = .IOI2 = 0.625, 
as = .OII2 = 0.375, aj = .III2 = 0.875, ag = .OOOI2 = 0.0625, etc. Note that the sequence 
alternates between the two half-unit intervals (0, 0.5) and (0.5, 1), which helps minimize 
spurious shock propagation. Furthermore the series can be shown to be optimal in terms 
of uniform coverage of the unit interval (0, 1). The value adopted in a given time step 
is the same for all inter-grid points. The exact solution at a given grid point is alternately 
taken to be either the left or right solution between two adjacent grid points. The solution 
is evaluated in each alternating half-time step at a time (1/2) At, where At = ucfl/ Ax 
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with CFL (Courant, Priedrichs, & Lewy 1967) number ticfl = 0.5. Thus the pure Ghmm's 
Method effectively adopts a CFL number of 0.5 for the full time step. Although most of 
our results use a simple one dimensional Cartesian grid, Wen et al. (1997) also present 
geometrical correction terms for carrying out one dimensional calculations in cylindrical or 
spherical symmetry. 

An important advance since Wen et al. (1997) are the studies generalizing the relativistic 
Riemann solution to include tangential flow (Pons, Marti, & Miiller 2000, Rezzolla, Zanotti, 
& Pons 2003). This allows one to extend Glimm's method to problems involving shear, and 
to begin to envision a two dimensional Glimm's method. Pons et al. obtain a solution by 
solving (1) the jump conditions across shocks, and (2) a differential equation that comes 
from a self-similarity condition along rarefaction waves. Rezzolla et al. present an integral 
solution to the equation derived by Pons et al. and they propose an efficient Gaussian 
quadrature technique for solving it. To solve the local Riemann problem between adjacent 
grid points we use the publicly available code RIEMANN_VT.F written by J.-M. Marti and E. 
Miiller (cf. Marti & Miiller 2003) which uses the formalism described in Pons et al. (2000) 
and Rezzolla et al. (2003). 



4. Testing 

Shock tube problems used in testing hydrodynamical codes are a subset of the Riemann 
problems, for which v = for all x. One dimensional Riemann problems are typically run on 
a grid such that < a; < 1, and the thermodynamic variables P, p, and v are discontinuous 
across x = 0.5 initially. Starting the simulation is equivalent to removing a diaphragm 
between left (L) and right (R) states. The strong gradients across x = 0.5 result in four 
constant states separated by three elementary waves: rarefaction, contact discontinuity, 
and shock wave. Analytical solutions for the time evolution of these problems for (special) 
relativistic hydrodynamics are given by Marti & Miiller (1994) for nonshearing problems, 
and by Pons et al. (2000) for Riemann problems with added shear (i.e., non-zero v±). 

The level of agreement between the exact, analytical solutions and the numerical ones 
is quantified by the Li norm error, defined for ID problems as Li = T,jAxj\uj — u{xj)\, 
where Xj is the coordinate of grid point j, u{xj) is the analytical value, and Uj the numerical 
value. The grid spacing is A.Xj. For consistency with previous groups, we take the solution 
in proper density. The analytical and numerical solutions are calculated on the same grids, 
and the number of grid points N in the solutions are varied between trials. 
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4.1. Riemann Problem 1 

The values in the initial left and right states are {p, p, v)l = (40/3, 10, 0) and 
(p, p, v)r = ((2/3) X 10-^ 1, 0). The adiabatic index T = 5/3. The result at t = 0.4 
is compared to the analytical one. The gradient in pressure p produces in the subsequent 
evolution a rarefaction wave moving left and a shock wave moving right, with a contact 
discontinuity between. The flow is mildly relativistic, with post-shock velocity v = 0.714. 
Figure [2] shows a comparison of the Glimm solution with the exact one, computed on a grid 
with N = 400. The small inset panels show a detail of the leading and trailing edges of the 
density spike associated with the shock. For the time step shown, the leading edge of the 
Glimm solution is off the analytical solution by one grid point, and the trailing edge is exact. 
Table [1] presents the Li errors in density between 3 methods, FLASH (from Morsony et al. 
2007), WENO (weighted essentially nonoscillatory, from Zhang & MacFadyen 2006), and 
Glimm. (This test problem has been studied by many workers previously — e.g., Hawley, 
Smarr, & Wilson 1984, Schneider et al. 1993, Marti & Miiller 1996, Wen, Panaitescu, & 
Laguna 1997, Marti et al. 1997, Aloy et al. 1999.) The asterisked values indicate those 
trials for which the leading and trailing shock edge positions of the Glimm solutions are in 
agreement with the analytical ones. 

For a small sample of individual Glimm trials, the Li error is not always a consistent 
indicator of success. Although shock propagation speeds are expected to be accurate in 
an averaged sense, within a given time step specific features in the Glimm solution can be 
one or two grid points off from their correct location. For problems with sharp edges, such 
as shocks, the error will be large (locally) at such a position. Most the rest of the error 
is introduced by idealizing the curved state (Riemann fan) to be composed of a series of 
piecewise constant states. Even if a shock edge location is incorrect at a given time step, 
at a slightly later time step, or at the same time step for a run with a different number of 
grid points A^, the Glimm solution may have the correct location of the shock front edges. 
Therefore a better way to measure the success of the method is to plot the Li errors for 
a large number of different trials, all compared at the same time step with the analytical 
solution for the same A^. For problems which are typically dominated by one large density 
enhancement, one observes bands of solutions representing those for which the calculated 
edges are (1) exact, (2) off by one grid point (leading or trailing edge), (3) off by two grid 
points total, (4) off by three grid points total, etc. We denote the cumulative grid point 
error in shock front localization by s. 

This effect is shown in Figure |3] where we plot the Li error versus A^. The black and 
blue points indicate values for the 6 grid points shown in Table [1], and the red points show a 
much larger sample drawn from ~ 10^ equi-logarithmically spaced A^ values for the Glimm 
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solutions. There is a large scatter vertically in the Glimm Li errors according to the degree 
of matching of the shock edges. 

4.2. Riemann Problem 2 

Riemann problem 2 has a more extreme pressure contrast between the L and R states 
initially than problem 1, and therefore drives a faster and higher density shock. The values 
in the initial left and right states are (p, p, v)l = (10^, 1, 0) and (p, p, v)r = (0.01, 1, 0). 
The adiabatic index F = 5/3. The result at t = 0.4 is compared to the analytical one. The 
flow is relativistic, with post-shock velocity v = 0.96. The shock speed is 0.986. The width 
of the shock is 6xs ~ 0.01 at t = 0.4 and for = 400 is covered by 4.2 grid points (in the 
analytical solution). The asterisked values indicate those trials for which s = 0. Table 2 
compares the Li errors for the three methods. 

Figure H] shows the Li error plot for Riemann problem 2, with the values given in 
Table [2] plus Glimm values for ~ 10^ additional N values. Due to the thinness of the shock 
compared to problem 1, there is now a clear banded structure to the Glimm solutions. The 
lowest striation, which also contains the first and sixth values from Table [2l corresponds to 
solutions for which both leading and trailing shock edge positions are exact, s = 0. The next 
highest striation, containing Glimm entries 3 — 5 from the table, corresponds to s = 1, and 
the third striation, containing the second Glimm entry from the table, corresponds to s = 2. 
The first striation lies about two orders of magnitude below the F errors, while the second 
and third are within a factor ~ 3 — 10 of F. 



4.3. Riemann Problem 3 

Riemann problem 3 starts with a strong negative pressure gradient that launches a 
reverse shock, and a positive flow speed in the left state that initiates a forward shock. Thus 
there is no Riemann fan. The values in the initial left and right states are (p, p, v)l = 
(1, 1, 0.9) and {p, p, v)fi = (10, 1, 0). The adiabatic index F = 4/3. The result at t = 0.4 
is compared to the analytical one. Table [3] compares the Li errors for the three methods. 

Figure [5] shows the Li error plot for Riemann problem 3, with the values given in Table 
|3]plus Glimm values for ~ 10^ more N values. None of the values in the table lie in the band 
for s = 1. The three upper limit triangles indicate solutions for which s = 0, i.e., all three 
shock edge locations in the problem are exact at t = 0.4. The only limiting precision is the 
machine epsilon e (~ 10~^^). For one dimensional problems consisting only of constant states. 
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Glimm's method finds the exact values, therefore the only error is introduced by shock edge 
location inaccuracies. In traditional methods this test problem produces postshock pressure 
oscillations in the reverse shock (e.g., Lucas-Serrano et al. 2004, see their Fig. 1; Zhang 
& MacFadyen 2006, see their Fig. 3). Lucas-Serrano et al. (2004) note, however, that the 
oscillations completely disappear when the CFL number is reduced below 0.3. 

4.4. "Easy" Shear: Riemann Problem 2 with {v±)r ^ 

We now proceed to one dimensional problems involving shear. The "easy" shear problem 
takes Riemann problem 2 and adds constant background shear in the R state, (f_L)_R = 0.99. 
The adiabatic index F = 5/3, and the result at t = 0.4 is compared to the analytical one. 
The highest Lorentz factor in the resulting flow 7 ~ 7.L Unlike purely Newtonian flows 
in which orthogonal components of the velocity field are decoupled from each other (aside 
from dissipation), with special relativity we now add the condition that ^ v\ < 1. This 
effectively limits the component of velocity along the direction of the flow v, and also the 
degree of density enhancement relative to background within the shock. In addition, 7 now 
includes a contribution from the shear. There is also a back reaction in terms of the evolution 
of v{x,t) on the initially constant v± values. Table H] compares the Li errors for the three 
methods. 

Figure [6] shows the Li errors for the values given in Table H plus ~ 10^ additional 
values for G. As with Figs. H] and [5l the banded structure associated with the precision in 
the shock edge localization is evident. The locus of solutions for s = lies ~ 10^ — 10'^ below 
the F errors, while the second striation, corresponding to s = 1, lies within a factor of 10 of 
the F errors. 

4.5. "Hard" Shear: Riemann Problem 2 with {v_i_)R 7^ and {v_i_)l ^ 

The "hard" shear problem starts with Riemann problem 2 and adds background shear in 
both the R and L states, {v±)r = (w_l)l = 0.9. The adiabatic index F = 5/3, and the result 
at t = 0.6 is compared to the analytical one. The highest Lorentz factor in the resulting 
flow is 7 ~ 35.8. Table [5] compares the Li errors for the three methods. The asterisked 
value indicates the trial for which s = 0. This problem poses a severe challenge for the 
traditional methods, but is well-handled by the Glimm method. In fact, the Li error for F 
for the highest values shown are equal to those for the lowest A^ values for G. Zhang & 
MacFadyen (2006) present results of the hard shear test for up to 51,200 grid points, either 
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uniform or the adaptive mesh equivalent (see their Table 7 and Fig. 9). Their Li errors for 
= 51200 of ~ 10~^ are comparable to those in our test for = 400. The challenge of 
relativistic ID shearing problems for standard FD techniques is also evident in Morsony et 
al. (2007, see their Fig. 24). The profiles of p and v for the FD shearing experiments shown 
in Mignone, Plewa, & Bodo (2005), Zhang & MacFadyen (2006), and Morsony et al. (2007) 
all exhibit a strong displacement and skewing of the shock density spike with respect to the 
analytical solutions. 

Figure [7] shows the Li errors for the values given in Table El plus ~ 10^ additional 
values for G. The s = striation lies ~ 10^ — 10^ below the F errors, and the higher striations 
are still a factor ~ 10 below F. 



4.6. Isentropic Smooth Flow 

4-6.1. Continuous Isentropic 

The previous problems contained sharp gradients produced by shocks. We now look at a 
problem with smooth flow, the isentropic flow problem. This consists of an initial state with 
smooth profiles in p, p, and v. A pulse of moving fluid is superposed on top of a constant 
density, zero velocity state. The velocity of each individual element is constant in time. 
Therefore the "exact" solution at a later time t > is found by advancing each element in 
time at its known velocity, which yields a grid with irregular spacing, and then interpolating 
the result back onto a uniform grid. 

The initial structure is given by 

Po{x)=p*[l + af{x)], (1) 

where p* is the density of the constant background state, and the function f{x) = (x^L~^ — 1)^ 
for |x| < L, and f{x) = for |x| > L. The width of the pulse is L and the amplitude is 
a. The initial velocity profile within the pulse is set by taking one of the two Riemann 
invariants to be constant, 

where cj. = Tp/{p + [r/(r — l)]p). The other Riemann invariant is not constant, 

. 1 ,„ (l±^) + 1 ,„ + . (3) 
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One inverts the equation for J_ to find the velocity 



(4) 



V = 



e^9 + V 



where 



9 = J- + 



1 



In 



) 



(5) 



Following previous workers (Zhang & MacFadyen 2006, Morsony et al. 2007) we use a 
domain —0.35 < x < 1, and adopt p* = 100, p* = 1, and v* = 0. We also take a = 1 
and L = 0.35. The adiabatic index F = 5/3, and the result at t = 0.8 is compared to the 
analytical one. Figure [8] shows the evolution of p, p, and v from the initial state. Table [6] 
compares the Li errors for the three methods. Figure [9] shows the Li errors for the values 
given in Table [6|, plus ~ 10^ additional N values for G. 



The Riemann problem and isentropic fiow problem span extremes of two possible initial 
states, one with constant states and one with smooth fiow. A better metric for realistic 
problems, where discontinuities and smooth fiows are found together, would combine these. 
Therefore we investigate the evolution of a structure that is initially piecewise isentropic: 
between the two isentropic parts we introduce a discontinuous jump in pressure and velocity. 
Since there is now no analytical solution, we carry out one ultra-high resolution run as the 
reference solution. 

The one change we make to the isentropic fiow problem is to force a jump in p at a; = 
such that the excess above the fioor level p = 100 drops by a factor of two. The sharp 
negative gradient in p at a; = drives a strong fiow to the right which is superposed on the 
natural fiow. Figure [TO] shows the evolution of p, p, and v from the initial state, and Figure 
[TT] shows the associated errors. The "exact" solution is obtained by computing a Glimm 
run for = 10^, and then interpolating to the grid spacing of each of the ~ 10^ trial runs. 
Since this is a modification of a standard test, there are no FD model errors with which to 
compare. 



In their generalization of the exact special relativistic Riemann problem to include shear. 
Pons et al. (2000) introduce a suite of 9 tests involving shear, also based on Riemann problem 



4-6.2. Piecewise Isentropic 



4.7. Shear Suite of Problems from Pons et al (2000) 
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2. These have been examined by Mignone, Plewa, & Bodo (2005) using the FLASH code 
(see their Fig. 5). In Figure [12] we present the resuhs of applying Ghmm's method to this 
test suite. As with the non-shearing test problems, constant states are reproduced exactly 
(i.e., to within machine precision), thereby avoiding the problems with FD methods alluded 
to earlier. 

4.8. Ultrarelativistic Shear Problems from Aloy & Rezzolla (2006) 

Rezzolla, Zanotti, & Pons (2003) study the effect of shear on the standard Riemann 
problems, and find that the standard pattern of a contact discontinuity sandwiched be- 
tween a rightward moving forward shock and a leftward moving reverse shock, abbreviated 
^SCS^, can be fundamentally altered by the presence of a strong shearing field. For suffi- 
ciently large shear, the reverse shock can be replaced by a rarefaction wave, hence the new 
pattern ^RCS^ arises. Aloy & Rezzolla (2006) explore the astrophysical ramifications of 
the Rezzolla et al finding as a potential mechanism for accelerating jets from AGNs, micro- 
quasars, and GRBs to very high Lorentz factors. They show that by varying the left hand 
pressure in a Riemann problem, one can change the nature of the solution. 

We present two additional shearing tests that delve deeper into the ultrarelativis- 
tic regime than the "hard" shear problem presented earlier. For the first case we take 
{p, p, V, 7)l = (10-3, 10-^ 0.99, 20) and (p, p, v, 7)/? = (10"^ lO'^, 0, 1). The Lorentz 
factor 7 includes both the normal and perpendicular velocities 7 = (1 — f ^ — t;^)^^/^. The 
adiabatic index F = 4/3, corresponding to the ultrarelativistic case. For this trial the shock 
speed Vs = 0.151. The result at t = 1.8 is compared to the analytical one. According to Aloy 
& Rezzolla (see their Fig. 4), = 10^^ should lie below the transition point from ^SCS^ 
to ^RCS^. Figure [13] shows a comparison between the Glimm's Method solution and the 
exact solution for N = 400, and figure UM shows the LI norm density errors at t = 1.8. Since 
this problem is relatively new, there are no published FD results with which to compare, but 
one suspects that the FD errors would be comparable or worse to those shown previously in 
connection with the "hard" shearing problem. 

For the second Aloy & Rezzolla shear problem we increase pi by eight orders of magni- 
tude to 10^. All other initial L and R parameters are the same. This p^ value should shift 
the wave pattern for the Riemann solution well into the regime ^RCS^ and yield a flow 
with maximum 7 ^ 10^ (Aloy & Rezzolla 2006 — see their Fig. 4). For this trial the shock 
speed Vs = 0.200. The result at t = 0.8 is compared to the analytical one. Figure [15] shows 
a comparison between the Glimm's Method solution and the exact solution for = 400, 
and figure [TB] shows the LI norm density errors at t = 0.8. For large N the Glimm solutions 
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acquire a permanent offset error in shock edge localization, rather than deviating about a 
mean s = 0. As with Fig. [H] we have only Glimm errors to present because the test is too 
new to have undergone published FD testing. 

4.9. Spherical Blast Wave 

The evolution of a relativistic blast wave in spherical symmetry has been examined by 
many workers. Panaitescu et al. (1997) present a detailed study using a hybrid Glimm/FD 
code, and taking 70 = 10^. Kobayashi & Zhang (2007) utilize a spherically symmetric 
relativistic code which uses a second-order Godunov method with an exact Riemann solver 
(described in Kobayashi, Piran, & Sari 1999) to investigate the evolution of a relativistic blast 
wave. Kobayashi & Zhang investigate a thin-shell case taking 70 = 10^, and a thick-shell 
case taking 70 = 10^. 

The final test shown in Wen et al. (1997) is that for a relativistic blast wave with 
initial Lorentz factor 70 = 10. For comparison in Figure [17] we show results for a run 
with similar starting conditions. To adapt to spherical geometry we use the geometrical 
correction terms given in Wen et al. (1997) with a = 2. Within a narrow radial range O.Olro 
centered at tq we initialize using a Blandford-McKee profile po{r) = 10^7qX~''^^7~^, where 
X = 1 + 16(1 — r/ro)7o, 7 ~ loX~^^'^^ 7o = 15, and tq = 0.4. We takepo(^) = 0.2pQ{r). Inside 
the initial shell po = Po = 10""^; outside the initial shell po = 1 and po = 10^"^. 

The profiles shown in Kobayashi & Zhang (2007) do not display obvious oscillations 
in the shocked shell. In our case, using a much smaller initial Lorentz factor, we see in 
Figure [T7I a number of small oscillations, particularly in 7. This indicates that the treatment 
of spherical geometry is worse than that of FD conservative methods such as the one of 
Kobayashi & Zhang. In addition, due to the sharpness of the density shell and the strong 
mass jumps accompanying grid points entering into and then leaving the shell, mass is 
conserved for the run shown in Figure [T71 only to within ~10%. 

5. Discussion 

We have presented the results of a series of tests done on standard problems in relativistic 
hydrodynamics using Glimm's method. To compare to previous works we utilize the Li 
norm errors in density. For problems involving smooth gradients such as the isentropic flow 
problem, Glimm's method fares worse than the standard finite difference techniques, due to 
the fact that solutions are typically off by ~ 1 — 2 grid points. In one dimension, however. 
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the constant states are exact to within machine precision. This is true irrespective of the 
presence of shear, thereby giving the method an advantage over FD methods. If there were 
only constant states in a solution, and if the leading and trailing shock edge locations were 
correct, then the entire solution would also be correct (to within machine precision). The 
idealization of piecewise constant states for the Riemann fan, however, is a source of error, 
as is the incorrect position of a shock edge. A better visualization of the Glimm errors than 
a simple table of Li errors versus grid point number N is achieved by calculating a large 
number of numerical and analytical values for varying N, and plotting the results. In such 
a plot one sees several bands of solutions corresponding to the total number of grid points 
s by which the shock edge locations are off. For a given problem, the degree to which sharp 
edges differ from their correct locations varies both with time within a given trial, and with 
A^. Therefore one cannot choose a priori the "right" resolution for any problem such that the 
errors arc minimized; one can only see what the errors are for being off the correct solution 
by a given s value. 

For the specific problems studied in this work, Riemann problem 1 yields similar global 
errors between Glimm and FD methods for the ensemble of ~ 10^ solutions. For Riemann 
problem 2, the Glimm errors are comparable to FD for solutions for which s ~ 3 — 4. The 
solutions with zero localization error s = (i.e., exact matching of the shock edges to their 
correct values) have Li errors ~ 10^ times smaller than the FD methods. For Riemann 
problem 3, the s — solutions are limited only by the machine e error, solutions for which 
s = 1 lie a factor ~ 10 below FD, and solutions with s ~ 2 — 4 are comparable to FD. For the 
easy shear problem, the s = solutions have errors ~ 10^ times smaller than for FD. The 
errors become comparable for s ~ 3 — 4. For the hard shear problem, the s — solutions 
have errors ~ 10^ — 10^ times smaller than for FD. The errors do not become comparable 
for any s. In fact, the Glimm errors for the lowest N values studied arc comparable to those 
for the highest N values in previous FD investigations. For smooth iscntropic flow, the FD 
errors are comparable to Glimm for the smallest values. For the largest A^ values, the 
FLASH errors are a factor ~ 10^-^ smaller than for Glimm, and for WENO ~ 10^-^ times 
smaller than Ghmm. For the relativistic blast wave test in spherical geometry (ID), the 
profiles are similar to those of a comparable run in Wen et al. (1997, see their Fig. 5). 

For the local Riemann problem, the Riemann solver RIEMANN_VT.F decomposes each 
solution into a left wave and a right wave. Depending on the conditions, many iterations 
may be required, therefore the computation time can varying greatly. Wen et al. (1997) 
discuss the slowness inherent in the Ghmm's Method and quote run times > 10 times slower 
than standard FD methods. Wc find, using a ~ 2GHz machine that, for example, Riemann 
problem 1 for A^ = 400 and t = 0.4 (640 half time steps) requires 7s of CPU time (27/is 
per grid point per half-time step), the hard shear problem for A" — 400 and t — 0.6 (960 



- 13 - 



half time steps) uses 17s (44/is per grid point per half-time step), and the Aloy & Rezzolla 
problem 2 for = 400 and t = 0.8 (1280 half time steps) takes 230s (450/is per grid point 
per half-time step). The = 10^ piecewise isentropic run required 4 wks. 

Although Glimm's method is superior in resolving shocks, for problems containing thin 
features, as is common in relativistic hydrodynamics, there is still a strong need for adaptive 
mesh refinement. For a given grid spacing, features are often too narrow to be resolved. 
Figure [TS] shows the variation of total grid mass m (computed from the proper density) with 
time for Riemann problem 2 for the six Glimm runs indicated in Table [21 (The ending time 
t = 0.4 is that for which the errors were calculated.) The abrupt vertical excursions in m arise 
as the shock widens with time and new grid points are incorporated into the shock feature. 
Since the density is higher within the shock, the mass jumps. For the higher A^ values there 
are always enough grid points to cover the shock, and the variation in total mass is small as 
the new shock grid points come into existence. For the lower A^ values, however, this is not 
the case. In fact, for the A^ = 100 run, there are no grid points representing the shock feature 
until t ~ 0.3, at which time the shock has widened to of order the grid spacing, and one grid 
point appears at the shock location, hence the large jump in mass. For Glimm's method to 
be a useful research tool, it will probably be necessary not only to have a two dimensional 
version, but also to include a provision for adaptive mesh refinement. Preliminary work 
on a 2D version has been encouraging, but more effort is required to address the issue of 
numerical stability. 



6. Conclusions 

We present the results of relativistic hydrodynamical tests using Glimm's method, along 
with a comparison to results using standard methods. Glimm's method in one dimension 
is superior to standard finite differencing for problems containing shocks, in which a sharp 
gradient appears. The introduction of shear does not degrade the quality of the solutions. 
Indeed, the work of Pons et al. (2000) generalizing the relativistic Riemann solution to 
include shear now also provides impetus for making a two dimensional relativistic Glimm's 
method. For problems involving smooth flow, the standard finite differencing methods are 
much better. Although constant states are calculated exactly (i.e., to within machine pre- 
cision) in Glimm's method, curved states such as Riemann fans are somewhat imprecisely 
modeled as being composed of a sum of piecewise constant states. Furthermore, the fact 
that there is an uncertainty of 1 — 2 grid points in the location of a given feature means 
that for models with smoothly varying physical parameters, the entire profile can be shifted 
slightly, leading to large global errors in comparison to an exact solution. The results of the 
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piecewise isentropic run indicate that for realistic applications containing both smooth flows 
and sharp gradients, standard FD methods give superior global behavior. Glimm's method 
may prove better for applications such as GRB afterglow shock propagation into a uniform 
medium where one is primarily interested in the physical evolution of high entropy material 
only within a restricted volume (i.e., the shocked gas), and not the global evolution of low 
density, low entropy regions far away from the shock. 

We thank Alin Panaitescu for helpful discussions concerning Glimm's method, and for 
allowing us to use the driver code from Wen, Panaitescu, & Laguna (1997) that sets up the 
hydrodynamical model and calls the Riemann solver. As mentioned earlier, we use the code 
RIEMANN_VT.F written by Jose Marti and Ewald Miiller for solving the Riemann problem. 
We also thank Brian Morsony for useful advice and a short IDL code to advance the exact 
solution for the isentropic flow problem, and Tod Strohmayer for a useful suggestion. Thanks 
also go to the anonymous referee for suggesting the piecewise isentropic flow problem and 
the ultrarelativistic shearing problems from Aloy & RezzoUa (2006). 
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Table 1. LI Error — Riemann Problem 1 



N 


FLASH 


WENO 


Glimm 


100 


0.13 


0.13 


0.029* 


200 


0.070 


0.074 


0.034 


400 


0.036 


0.033 


0.017 


800 


0.018 


0.021 


0.0035* 


1600 


0.0085 


0.010 


0.0033 


3200 


0.0043 


0.0051 


0.0069 
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Table 2. LI Error — Riemann Problem 2 



N 


FLASH 


WENO 


Glimm 


100 


0.21 


0.21 


0.0034* 


200 


0.15 


0.14 


0.10 


400 


0.083 


0.093 


0.024 


800 


0.046 


0.055 


0.012 


1600 


0.025 


0.025 


0.0061 


3200 


0.013 


0.015 


0.00011* 
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Table 3. LI Error — Riemann Problem 3 



N 


FLASH 


WENO 


Glimm 


100 


0.059 


0.10 


0.061 


200 


0.035 


0.063 


0.031 


400 


0.021 


0.030 


0.013 


800 


0.013 


0.017 


0.0070 


1600 


0.085 


0.095 


0.0038 


3200 


0.033 


0.052 


0.0019 
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Table 4. LI Error — Easy Shear 



N 


FLASH 


WENO 


Glimm 


100 


0.63 


0.76 


0.24 


200 


0.34 


0.39 


0.12 


400 


0.17 


0.23 


0.059 


800 


0.084 


0.12 


0.029 


1600 


0.044 


0.066 


0.015 


3200 


0.023 


0.034 


0.029 
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Table 5. LI Error - Hard Shear 



N 


FLASH 


WENO 


Glimm 


100 


0.51 




0.038 


200 


0.46 




0.019 


400 


0.33 


0.52 


0.0096 


800 


0.22 


0.36 


0.00048* 


1600 


0.13 


0.23 


0.0030 


3200 


0.083 


0.13 


0.0029 


6400 


0.053 


0.065 


0.0013 
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Table 6. LI Error — Isentropic Flow 



N 


FLASH 


WENO 


Glimm 


80 


5.5e-3 


2.1e-3 


0.0072 


160 


1.6e-3 


l.lc-4 


0.0052 


320 


4.0e-4 


1.7C-5 


0.0033 


640 


l.Oe-4 


1.5e-6 


0.0024 


1280 


2.5e-5 


1.6e-7 


0.0019 


2560 


5.4e-6 


1.9^8 


0.0014 


5120 


1.6e-6 


2.4e-9 


0.00053 
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Fig. 1. — The evolution of Lorentz factor 7 for a Blandford-McKee initial state with 70 = 5 
in a 3D Cartesian calculation using the method described in del Zanna & Bucciantini (2002), 
using the local Lax-Friedrichs flux. The four panels show increasing grid resolution in a slice 
along the propagation direction. The initial step is the leftmost profile in each panel, and 
profiles moving to the right show the shock development at eight subsequent time steps. For 
ease of viewing, the solutions in the first panel are connected by sohd lines. The dotted 
curve in each panel indicates the 7 value corresponding to the local maxima in p for the nine 
time steps. (For the first two panels there is a [spurious] offset between the local maxima in 
p and 7.) The number of grid points per small tick mark is {top to bottom) 1, 4, 16, and 64, 
respectively. 

Fig. 2. — A comparison of the pressure p {open triangles), density p {open squares), and 
velocity v {open pentagons) for Riemann problem 1 at t = 0.4 using Glimm's method with 
the exact solution {solid line), where both are computed on a grid with = 400. The small 
panels show the detail of the leading and trailing edges of the density spike. A dashed line 
connects the Glimm density points. 

Fig. 3. — The Li errors in density for the F (FLASH and WENO — shown in black) and G 
solutions (Glimm) in Table [1]— Riemann problem 1 (shown in blue), as well as the results 
using ~ 10^ additional A^ values for G (shown in red). 

Fig. 4. — The Li errors in density for the F (FLASH and WENO — black) and G (Glimm 

— blue) solutions in Table [2] (Riemann problem 2), as well as ~ 10^ additional A^ values for 
G (red). The fact that the shock is narrower than for Riemann problem 1 leads to a more 
pronounced striationing; with fewer points spanning the shock, the relative error introduced 
by being off a given number of grid points in the shock edge location is larger. 

Fig. 5. — The Li errors in density for the F (FLASH and WENO — black) and G (Glimm 

— blue) solutions in Table [3] (Riemann problem 3), as well as ~ 10^ more A^ values for G 
(red). The solutions listed in the table all lie in the band for which s = 2. The upper limit 
triangles indicate three solutions which are limited only by machine precision (~ 10^^^). 

Fig. 6. — The Li errors in density for the F (FLASH and WENO — black) and G (Glimm 

— blue) solutions in Table H] ("easy" shear), as well as ~ 10^ more A^ values for G (red). 

Fig. 7. — The Li errors in density for the F (FLASH and WENO — black) and G (Glimm 

— blue) solutions in Table [5] ("hard" shear), as well as ~ 10^ more A^ values for G (red). 

Fig. 8. — The evolution of p {top panel), p {middle panel), and v {bottom panel) for the isen- 
tropic flow problem, taking A^ = 400. Shown are the initial state {dotted) and 5 subsequent 
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equally spaced times steps up to t = 0.8. 

Fig. 9.— The Li errors in density for the F (FLASH - black), W (WENO - black), and G 
(Ghmm — blue) solutions in Table 6 (isentropic flow), as well as ~ 10^ additional N values 
for G (red). 

Fig. 10. — The evolution of p {top panel), p {middle panel), and v {bottom panel) for the 
picccwise isentropic flow problem, taking = 400. Shown are the initial state {dotted) and 

5 subsequent equally spaced times steps up to t = 0.8. 

Fig. 11. — The Li errors in density for the Glimm solutions for the piecewise isentropic flow 
problem. 

Fig. 12. — Comparison of numerical with analytical solutions for the shearing test problems 
shown in Pons et al. (2000, see their Fig. 4) and Mignone et al. (2005, see their Fig. 5) 
using N — 400. The small insets in each panel show details of the leading and trailing edges 
of the density spike. The tests begin with Ricmann problem 2, and then add, for the initial 
shear in the R and L states {from left to right), {v±)r = 0, 0.9, and 0.99, and {from top 
to bottom) {vj_)L = 0, 0.9, and 0.99. The polytropic index F = 5/3, and the solution is 
evaluated at t = 0.4. Thus the upper left panel shows the solution for Riemann problem 2 
from Section 4.2, the upper right panel shows the "easy" shear solution from Section 4.3, 
and the central panel shows the "hard" shear solution from Section 4.4 (evaluated eikt — 0.4, 
however, rather than t — 0.6). For the central and lower right panels, s = 0; for all other 
panels, s = 1. 

Fig. 13. — A comparison of the pressure p {open triangles), density p {open squares), velocity 

V {open pentagons), and Lorentz factor (= [1— w ^— w^]"^/^ — open hexagons) for the first Aloy 

6 RezzoUa ultrarelativistic problem at t — 1.8 with the exact solution {solid line), where 
both are computed with N — 400. The small insert panel shows a detail of the density spike. 
A dashed line connects the Glimm density points. 

Fig. 14. — The Li errors in density accompanying the Glimm solutions for the first Aloy & 
RezzoUa ultrarelativistic problem, at t = 1.8. 

Fig. 15. — A comparison of the pressure p {open triangles), density p {open squares), velocity 

V {open pentagons) , and Lorentz factor (7 = [1 — — f^]^^/^ — open hexagons) for the second 
Aloy & RezzoUa ultrarelativistic problem at t = 0.8 with the exact solution {solid line), where 
both arc computed with = 400. A snapshot of the first three variables is shown in the 
top panel, and the bottom panel shows 7. The small insert within the top panel presents a 
detail of the density spike. A dashed line connects the Glimm density points. 
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Fig. 16. — The Li errors in density accompanying the Ghmm solutions for the second Aloy 
& Rezzolla ultrarelativistic problem, at t = 0.8. 

Fig. 17. — The evolution of {from top to bottom) Lorentz factor 7, density p, and pressure 
p for a spherically symmetric ID test run of a thin-shell, relativistic blast wave to compare 
with Wen et al (1997, see their Fig. 5). For this run = 10^ over the entire computational 
domain (0.075 < r < 5.1), or 3800 grid points over the domain plotted. A Blandford-McKee 
profile with Lorentz factor 70 = 15 is taken initially for a thin spherical shell extending from 
0.99rs to Vg, where = 0.4. The frame of reference is continually adjusted so that the 
origin corresponds to the position of the contact discontinuity. There is a rightward moving 
forward shock and a leftward moving reverse shock. 

Fig. 18. — The variation of total mass with time, integrated over the grid, for Riemann 
problem 2. The six panels accompany the six Glimm entries in Table [2l 
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